% The purpose of this program is employed to construct Chebyshev polynomial approximation ...
% coupled with Smolyak sparse grid numerical integration method.The program can be verified through two examples from Reference 1.
% Author: W.Q. 2020.08.20
% Reference 1: Wu Jinglai. Dynamics Uncertainty Research Based on Interval Arithmetic Using Chebyshev Polynomials(Chinese), Phd thesis, 2013.[P111 ~ P119]
% Reference 2: Xiong Fenfen. Engineering Uncertainty Analysis by Probabilistic Method(Chinese),science press, 2015.[P98 ~ P104] 
tic;
syms x y                                                                        % Construction of symbolic variables: x and y
Polynomial_Order = 8;                                                           % The total degree of polynomials
Function_Handle = @Multi;                                                       % Function handle variable, such as 'Multi' and 'Schaffer'
Ndim = 2;
a=-1e-2;b=1e-2;det=1e-3;														% It is assumed that all variables have the same ranges; such as [a,b]
X = a:det:b;
Y = a:det:b;
[Coefficient_Set,Degree_Nd_Set] = Surrogate_Model_Nd(a,b,Polynomial_Order,Function_Handle);  % The coefficient array and degree distribution array corresponding to N dimensional Chebyshev Polynomials
Real_Result_Array = zeros(length(X)*length(Y),1);                               % Accurate function value array
Rough_Result_Array = zeros(length(X)*length(Y),1);                              % Approximate function value array
% Reduced_Rough_Result_Array = zeros(length(X)*length(Y),1);					% Simplified approximate function value array, which ignore some small coefficient elements
[ i_X,  i_Y ] = ndgrid(1:length(X),1:length(Y));								% N dimensional grid nodes
Index_Nd = [ reshape(i_X,[length(X)^Ndim,1]) reshape(i_Y,[length(Y)^Ndim,1])];

for index = 1 : size(Index_Nd,1)
		Rough_Result_Array(index) = Numerical_Surrogate_Model_Expression(...
			                       Coefficient_Set,Degree_Nd_Set,a,b,X(Index_Nd(index,1)),Y(Index_Nd(index,2)));
		Real_Result_Array(index) = Function_Handle(X(Index_Nd(index,1)),Y(Index_Nd(index,2)));
% 		Reduced_Rough_Result_Array(index) = ...
% 			- 41.329462452388554538629250600934*X(Index_Nd(index,1))^4 + 6.2831848481709196008182516379748*X(Index_Nd(index,1))^2 - ...
% 			41.329463273953592761245090514421*Y(Index_Nd(index,2))^4 + 6.2831848482552965506897635350469*Y(Index_Nd(index,2))^2 + ...
% 			1.00000000000510155806487944119;
end
Symbolic_Expression = Symbolic_Surrogate_Model_Expression(Coefficient_Set,Degree_Nd_Set,a,b,x,y);  % The symbolic expression of N dimensional Chebyshev Polynomials

% Result Visualization Section
[xx,yy] = meshgrid(X,Y);
figure(1);hold on;
plot3(xx,yy,reshape(Real_Result_Array,length(X),length(X)));					% 3D visualization of Accurate points
% plot3(xx,yy,reshape(Reduced_Rough_Result_Array,length(X),length(X)),'b*');
figure(2);hold on;
plot3(xx,yy,reshape(Rough_Result_Array,length(X),length(X)));					% 3D visualization of Approximate points

% Evaluation Section
Max_Error_Parameter = max(abs(Real_Result_Array-Rough_Result_Array))/...
                      max(abs(Real_Result_Array));								%  No. 1 Quantitative index from Reference 1
Mean_Error_Parameter = mean(abs(Real_Result_Array-Rough_Result_Array))/...
                       max(abs(Real_Result_Array));								%  No. 2 Quantitative index from Reference 1
toc;